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Abstract 



We present an object-oriented open-source framework for solving the dynamics of open quantum systems written in 
Python. Arbitrary Hamiltonians, including time-dependent systems, may be built up from operators and states defined 
by a quantum object class, and then passed on to a choice of master equation or Monte-Carlo solvers. We give an 
overview of the basic structure for the framework before detailing the numerical simulation of open system dynamics. 
Several examples are given to illustrate the build up to a complete calculation. Finally, we measure the performance 
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of our library against that of current implementations. The framework described here is particularly well-suited to the 
fields of quantum optics, superconducting circuit devices, nanomechanics, and trapped ions, while also being ideal for 
use in classroom instruction. 
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1. Introduction 

Every quantum system encountered in the real world is 
an open quantum system [T]. For although much care is 
taken experimentally to eliminate the unwanted influence 
of external interactions, there remains, if ever so slight, 
a coupling between the system of interest and the exter- 
nal world. In addition, any measurement performed on 
the system necessarily involves coupling to the measuring 
device, therefore introducing an additional source of ex- 
ternal influence. Consequently, developing the necessary 
tools, both theoretical and numerical, to account for the 
interactions between a system and its environment is an 
essential step in understanding the dynamics of quantum 
systems. 

By deflnition, an open quantum system is coupled to 
an environment, also called a reservoir or bath, where 
the complexity of the environmental dynamics renders the 
combined evolution of system plus reservoir intractable. 
However, for a system weakly coupled to its surroundings, 
there is a clear distinction between the system and its en- 
vironment, allowing for the dynamics of the environment 
to be traced over, resulting in a reduced density matrix 
describing the system alone. The most general dynamical 
equation governing this reduced system density matrix is 
given by the Lindblad master equation S2; describing the 
evolution of an ensemble average of a large (formally in- 
finite) number of identical system realizations. Although 
the density operator formalism sufficed for the first half- 
century of quantum mechanics, the advent of single-ion 
traps in the 1980's ^ motivated the study of the quantum 
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trajectories, or Monte-Carlo, description for the evolution 
of a single realization of a dissipative quantum system [3] . 

In general, for all but the most basic of Hamiltonians, an 
analytical description of the system dynamics is not pos- 
sible, and one must resort to numerical simulations of the 
equations of motion. In absence of a quantum computer 
[S], these simulations must be carried out using classical 
computing techniques, where the exponentially increas- 
ing dimensionality of the underlying Hilbert space severely 
limits the size of system that can be efficiently simulated 
[6l [7]. However, in many fields such as quantum optics 
[HI [H] J trapped ions [3J [TO] , superconducting circuit devices 
[TTl [T21 [T5] . and most recently nanomechanical systems 
[HI [TOl [TOl [TT] . it is possible to design systems using a 
small number of effective oscillator and spin components, 
excited by a small number of quanta, that are amenable 
to classical simulation in a truncated Hilbert space. 

Of the freely available quantum evolution software pack- 
ages [TOJ[TOl[2n], the Quantum Optics Toolbox (qotoolbox) 
[TO] has by far been the most successful. Although orig- 
inally geared toward quantum optics, the qotoolbox has 
gained popularity in a variety of disciplines, driven in part 
by its rapid modeling capabilities and easy to read code 
syntax. Yet, at the same time, the qotoolbox has not been 
updated in nearly a decade, leaving researchers to rely on 
an outdated numerical platform. Moreover, while the code 
underlying the qotoolbox is open-sourced, it does rely on 
the proprietary Matlab [21j computing environment mak- 
ing it an impractical solution for many research groups, as 
well as for use as an educational tool inside the classroom. 

In this paper, we describe a fully open-source imple- 
mentation of a framework designed for simulating open 
quantum dynamics written in the Python programming 
language [22J called the Quantum Toolbox in Python or 
QuTiP [23]. This framework distinguishes itself from the 
other available software solutions by providing the follow- 
ing advantages: 

• Based entirely on open-source software. 

• Easy to read, rapid code development using the 
Python programming language. 

• Support for arbitrary, time dependent Hamiltonians. 

• Makes use of the multiple processing cores found in 
modern computers. 

• Community based infrastructure, allowing for user 
contributions to the code base. 

Although Python is an interpreted programming lan- 
guage, it is well suited for scientific computations as a 
result of its large collection of high-performance low-level 
numerical libraries [24] . mathematical functions j25j, and 
data visualization capabilities [26) , that largely are imple- 
mented in efficient compiled code. In particular, QuTiP 



relies heavily on the sparse matrix and dense array func- 
tionality provided by the SciPy [2S] and NumPy [23] pack- 
ages, respectively. Since the bulk of a typical calculation 
is spent in these libraries, a QuTiP simulation can achieve 
nearly the same performance as compiled code. The ad- 
vantage of using the Python programming language over a 
compiled programming language is a greatly simplified de- 
velopment process, and more transparent, less bug-prone 
code. For data visualization QuTiP uses the matplotlib 
package [26], which is capable of producing publication- 
quality 2D and 3D figures in a wide range of styles and 
formats. 

Given the success of the qotoolbox, the development of 
QuTiP has in part been directed toward providing a re- 
placement for this excellent, yet aging software. In the 
spirit of open-source development, we have strived to use 
the best parts of the qotoolbox in QuTiP, while improving, 
replacing, or complementing the parts that were in need 
of modernization. The result is a framework for simulat- 
ing quantum system dynamics that is in many ways more 
efficient and better suited for modern computers, as well 
as better positioned for further development and adoption 
to new computer architecture advances. Given the size of 
the QuTiP framework, we do not hope to cover all of its 
functionality here. Instead, we will focus on the key data 
structures, and numerical routines underlying the major- 
ity of calculations. In addition, we will highlight a variety 
of example calculations that we hope will give the reader 
a flavor of the capabilities of QuTiP, and highlight what 
is possible using this framework. A complete overview of 
QuTiP is given on its website |23[ . 

This paper is organized as follows. In Sec. [2] we in- 
troduce the main QuTiP class, representing a quantum 
operator or state vector, and its associated data struc- 
tures and methods. In Sec. [3] we give a brief overview of 
the density matrix formalism before discussing the master 
equation and Monte-Carlo methods used in QuTiP. Sec- 
tion [4] presents a selection of examples meant to illustrate 
how calculations are performed using the QuTiP frame- 
work. Section [5] compares the performance of the QuTiP 
master equation and Monte-Carlo solvers to those in the 
qotoolbox. Finally, Sec. [6] briefly concludes, while a list 
of user accessible functions built into QuTiP, as well as 
example codes, are relegated to the appendix. 

2. The QuTiP framework 

QuTiP provides an object-oriented framework for rep- 
resenting generic quantum systems, and for performing 
calculations and simulations on such systems. In order 
to simulate a quantum system, one must first construct 
an object that encapsulates the properties of an arbitrary 
state vector or operator. A unified representation of quan- 
tum operators and state vectors is implemented in QuTiP 
by means of the quantum object class (Qobj), that uses 
a sparse matrix representation of a quantum object in a 
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Figure 1: (Color) A quantum object, operator or state vector, is 
represented by a quantum object class (Qobj) instance. The Qobj 
class may be thought of as a container, holding the data structures 
required to fully characterize a generic quantum object, as well as a 
list of instructions on how to manipulate these items. The primary 
data structures are the data, in sparse matrix form, that represents a 
quantum object in a given Hilbert space, the type of the object (ket, 
bra, operator, super-operator), whether it is Hermitian or not, and 
the objects dimension and shape. Here, the dimension describes the 
structure of the Hilbert space, i.e., whether it is a composite system 
and how it is composed. The Qobj class also defines a variety of 
methods that implement common functions operating on quantum 
objects. See Table [T] for a list of the Qobj class methods. 



finite dimensional Hilbert space. The Qobj class inter- 
nally maintains a record of the principal attributes of the 
quantum object it represents. These include, the objects 
type (i.e ket, bra, operator, or super-operator), whether 
the underlying object is Hermitian, the dimensionality of 
a composite object formed via the tensor-product, and the 
size of the sparse data matrix. A schematic illustration of 
the key components underlying the Qobj class is shown in 
Fig.g 

In addition to serving as a book-keeper for the prop- 
erties of a quantum object, the Qobj class is also a com- 
putational object, implementing the usual binary arith- 
metic operations, and a variety of class methods for per- 
forming common object manipulations as presented in Ta- 
ble [l] Therefore, with just a few lines of QuTiP code, it is 



Method 


Description 


dagO 


Adjoint of the quantum object. 


diagO 


Diagonal elements of object. 


eigenstatesO 


Eigenstates and eigenvectors. 


expm ( ) 


Exponentiated quantum object. 


fullO 


Dense array representation. 


normO 


L2 norm (states), trace norm (oper). 


sqrtmO 


Matrix square root. 


trO 


Trace of quantum object. 


unit 


Normalizes the quantum object. 



Table 1: List of methods built into the Qobj class. 



straightforward to construct Hamiltonians from arbitrary 
combinations of operators, and to construct density matri- 
ces and state vectors that represent complicated superpo- 
sitions of basis states. To further simplify this important 
step, QuTiP provides a library of commonly occurring op- 
erators and states which are given in [Appendix A[ 

For example, to create an instance of the Qobj class that 
represents the ubiquitous two-level Hamiltonian {h = 1) 

H = ]^e(T, + ^AfT,, (1) 

with energy splitting e and transition energy A, one can 
use the following QuTiP code: 

H = 0.5 * epsilon * sigmazO + 0.5 * delta * sigmaxO 

where epsilon and delta represent user defined con- 
stants. The result is a single Qobj instance H that rep- 
resents the Hamiltonian operator Eq. ([T]). 

Composite quantum systems are nearly as easy to cre- 
ate. Consider the Jaynes-Cumming Hamiltonian 

H — ojqo) a + 2^'^^ ^ di^'^ + o.'^-)^ (2) 

with cavity frequency ujq and coupling strength g, describ- 
ing a cavity field coupled to a two-level atom (qubit). 
A Qobj instance representing this composite Hamiltonian 
can be created with the following code: 

a = tensor (destroy(N) , qeye(2)) 
sm = tensor (qeye(N) , destroy(2)) 
sz = tensor (qeye (N) , sigmazO) 
H = omegaO * a.dagO * a + 0.5 * epsilon * sz 
+ g * (a.dagO * sm + a * sm.dagO) 

where the tensor function is used to construct compos- 
ite operators for the combined Hilbert space of the cavity 
(truncated to the N lowest Fock states) and the atom. 

Since the Qobj class provides a unified representation 
for operators and states, we can use exactly the same tech- 
nique to generate quantum states (either as state vectors 
or density matrices) . A possible initial state for the system 
described by Eq. (l2|, is generated with the QuTiP code: 

psiO = tensor(f ock(N,0) , (f ock(2,0)+f ock(2, D) .unitO) 

creating the Qobj representation of a cavity in its ground 
state, coupled to a qubit in a balanced superposition of its 
ground and excited state = (|0)^ |0), + |0)^ \1)^)/V2. 

Here, the subscripts c and q denote the cavity and qubit 
states, respectively. The normalization factor (\/2) is ap- 
plied automatically using the unitO method. 

In the previous example, we have used builtin QuTiP 
library functions to generate Qobj instances of commonly 
occurring operators and states, and the associated arith- 
metic operations in the Qobj clasfj^ to combine these 
quantum operators into more complicated systems. The 



■^The binary arithmetic operators +, -, and * are defined for two 
Qobj objects, and -|-, -, * as well as / are defined between a Qobj 
object and a real or complex number. 
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close correspondence between the mathematical formula- 
tion and the programming code makes it easy and trans- 
parent to define quantum objects. This is especially impor- 
tant when working with quantum systems of more complex 
structure. Note that the Qobj instances in the previous ex- 
amples are all self-contained descriptions of the quantum 
object they represent. From the Qobj instance alone, a 
number of properties pertaining to the quantum object 
may be calculated, and various operations be applied (see 
Table [I]). This includes, for example, operations such as 
the Hermitian adjoint, normalization, and trace, as well 
as computation of the eigenstates and eigenenergies and 
operator exponentiation. 

In addition, QuTiP includes several important func- 
tions opera ting on multiple states and/or operators (see 
Table A.2|. We have already seen one such example 
in the tensor function used to generate tensor product 
states. These states may also be decomposed into their 
constituent parts by performing a partial trace over se- 
lected degrees of freedom using the function ptrace. For 
example, from the composite wave function psiO for the 
oscillator-qubit system, the state of the qubit can be ex- 
tracted by tracing out the oscillator degrees of freedom 
using the QuTiP code: 

rhoO_qublt = ptrace (psiO, 1) 

where the second argument is the index of the system that 
we wish to keep. In general, it can be a list of indices. 
The properties of the resulting Qobj instance (shown in 
Fig. nj may be inspected by displaying the string repre- 
sentation of the object returned by its str method. This 
method is implicitly invoked when the object is printed to 
the standard output 



[2]] , shape = [2, 2] , 



print rhoO_qubit 
Quantum object: dims = [[2], 
type = oper, isHerm = True 
Qobj data = 
[[0.5 0.5] 
[0.5 0.5]] 



which, in this case, shows that the Qobj instance 
rhoO_qubit is a 2 x 2, Hermitian quantum operator rep- 
resenting a balanced coherent superposition of its two ba- 
sis states. From a Qobj instance one may also calculate 
items such as the expectation value (expect) for an arbi- 
trary operator with the QuTiP function, find the fidelity 
(fidelity) between two density matrices [27| . or calculate 
the Wigner function (wigner) of a quantum state. Using 



these, and other functions (Table. |A72), in the exploration 



of open quantum dynamics will be the focus of Sec. |4j 

Even though the emphasis of QuTiP is on dynamical 
modeling, it is possible to obtain nontrivial results directly 
from a quantum object. As an example, let us consider 
the Jaynes-Cummings model in the ultra-strong coupling 
regime g > ujQ,e where the rotating wave approximation 
(RWA) is no longer valid 



(3) 



0.0 




1.0 1.5 
Coupling strength ff/wo 



Figure 2: (Color) Expectation value for the number of excitations 
in the cavity (blue) and qubit (dashed-red) modes of the non-RWA 
Jaynes-Cummings model Eq. Jsl as the coupling strength g is in- 
creased into the ultra-strong coupling regime g/ojQ > 1. The inset 
figure displays the Wigner function for the cavity mode at the largest 
coupling strength, g = 2.5a;o, where ojq is the bare cavity frequency. 
At this coupling value, the state of the system is well-approximated 
by Eq. @. 



coupling strengths in superconducting circuit devices |31j . 
When the coupling strength g is a significant fraction of 
the cavity and qubit frequencies, the ground state of the 
cavity mode, after tracing out the qubit, is no longer the 
vacuum state. Instead, the anti-resonant terms propor- 
tional to a^cr_(. and ao- give rise to an anomalous ground 
state which, in the large coupling limit g/wg ^ 1, may be 
approximated as [311 US] 
1 

I q I I c\ I q 



1^9 



72 



(4) 



Recently, this regime has become of interest [IHl [2S1 130] 

due to the experimental realization of the required large 



where the cavity mode is in a Schrodinger cat-state with 
\a\ ~ g. This ground state can be evaluated by finding 
the eigenstates and eigenvalues of the Hamiltonian, and 
can therefore be extracted directly from the Qobj repre- 
sentation of Eq. ([3]). In Fig. ([2]) we plot the cavity and 
qubit occupation numbers for the groundstate of Eq. (|3]) 
as a function of the coupling strength. Here, the cav- 
ity is on resonance with the qubit transition frequency, 
Wo = e = 27r. In addition. Fig. [2] shows the Wigner func- 
tion for the cavity mode at the largest coupling strength 
g — 2.5a;o, which is well approximated by Eq. ([4|. The 20 
lines of QuTiP code used in calculating Fig. [2] are given in 
[Appendix B.l| 

3. Evolution of open quantum systems 

The main focus of QuTiP is the time-evolution of open 
quantum systems. Before we describe how this problem is 
approached in QuTiP, we give a brief review of the the- 
ory of quantum evolution, and the available methods for 
numerically integrating the equations of motion. 
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The dynamics of a closed (pure) quantum system is gov- 
erned by the Schrodinger equation 

where ^> is the wave function, H the Hamiltonian, and h 
is Planck's constant. In general, the Schrodinger equation 
is a partial differential equation (PDE) where both ^E* and 
H are functions of space and time. For computational 
purposes it is useful to expand the PDE in a set of basis 
functions that span the Hilbert space of the Hamiltonian, 
and to write the equation in matrix and vector form 

ihf^\^)=H\^), (6) 

where \^p) is the state vector and H is the matrix repre- 
sentation of the Hamiltonian. This matrix equation can, 
in principle, be solved by diagonalizing the Hamiltonian 
matrix H. In practice, however, it is difficult to perform 
this diagonalization unless the size of the Hilbert space 
(dimension of the matrix H) is small. Analytically, it is a 
formidable task to calculate the dynamics for systems with 
more than two states. If, in addition, we consider dissipa- 
tion due to the inevitable interaction with a surrounding 
environment, the computational complexity grows even 
larger, and we have to resort to numerical calculations in 
all realistic situations. This illustrates the importance of 
numerical calculations in describing the dynamics of open 
quantum systems, and the need for efficient and accessible 
tools for this task. 

While the evolution of the state vector in a closed quan- 
tum system is deterministic, open quantum systems are 
stochastic in nature. The effect of an environment on the 
system of interest is to induce stochastic transitions be- 
tween energy levels, and to introduce uncertainty in the 
phase difference between states of the system. The state 
of an open quantum system is therefore described in terms 
of ensemble averaged states using the density matrix for- 
malism. A density matrix p describes a probability distri- 
bution of quantum states \ijjn), in a matrix representation 
p = J2nPn IV'n) {ipnl, wherc p„ is the classical probability 
that the system is in the quantum state \tjjn)- The time 
evolution of a density matrix p is the topic of the remaining 
portions of this section. 

3.1. Master equation 

The standard approach for deriving the equations of mo- 
tion for a system interacting with its environment is to 
expand the scope of the system to include the environ- 
ment. The combined quantum system is then closed, and 
its evolution is governed by the von Neumann equation 

Ptot (^) = ^ ^ [^tot ' Aot (*)] ' C^) 

the equivalent of the Schrodinger equation ([s]) in the den- 
sity matrix formalism. Here, the total Hamiltonian 

-fftot — Hsys + Hcnv + ^^int , (8) 



includes the original system Hamiltonian -ffgysi the Hamil- 
tonian for the environment i?cnv, and a term represent- 
ing the interaction between the system and its environ- 
ment i?int- Since we are only interested in the dynamics 
of the system, we can at this point perform a partial trace 
over the environmental degrees of freedom in Eq. ([7|, and 
thereby obtain a master equation for the motion of the 
original system density matrix. The most general trace- 
preserving and completely positive form of this evolution 
is the Lindblad master equation for the reduced density 
matrix p = Tronv[ptot] 

-^^[Hit),pit)] 

+ E ^ [2C„p(t)C,t - p{t)CiCn - CiCnPit)] , (9) 

n 

where the C„ — ^/%An are collapse operators, and An 
are the operators through which the environment couples 
to the system in i?int, and 7„ are the corresponding rates. 
The derivation of Eq. ^ may be found in several sources 
[H [331 [Ml , and will not be reproduced here. Instead, we 
emphasize the approximations that are required to arrive 
at the master equation in the form of Eq. and hence 
perform a calculation in QuTiP: 

Separability: At i = there are no correlations be- 
tween the system and its environment such that the 
total density matrix can be written as a tensor prod- 
uctp(„t(0) = /(0)<8pLv(0). 

Born approximation: Requires: (1) that the state 
of the environment does not significantly change as a 
result of the interaction with the system; (2) The sys- 
tem and the environment remain separable through- 
out the evolution. These assumptions are justified 
if the interaction is weak, and if the environment is 
much larger than the system. In summary, ptot(^) ~ 

p{t) (g) Ponv 

Markov approximation: The time-scale of decay 
for the environment r^nv is much shorter than the 
smallest time-scale of the system dynamics Tgys ^> 
Tenv This approximation is often deemed a "short- 
memory environment" as it requires that environmen- 
tal correlation functions decay on a time-scale fast 
compared to those of the system. 

Secular approximation: Stipulates that elements 
in the master equation corresponding to transition fre- 
quencies satisfy jwab— cjcdl ^ l/'''sys: 1-6., all fast rotat- 
ing terms in the interaction picture can be neglected. 
It also ignores terms that lead to a small renormaliza- 
tion of the system energy levels. This approximation 
is not strictly necessary for all master-equation for- 
malisms (e.g., the Block-Redfield master equation), 
but it is required for arriving at the Lindblad form 
^ which is used in QuTiP. 



5 



For systems with environments satisfying the conditions 
outhned above, the Lindblad master equation ^ governs 
the time-evolution of the system density matrix, giving 
an ensemble average of the system dynamics. In order 
to ensure that these approximations are not violated, it 
is important that the decay rates 7„ be smaller than the 
minimum energy splitting in the system Hamiltonian. Sit- 
uations that demand special attention therefore include, 
for example, systems strongly coupled to their environ- 
ment, and systems with degenerate or nearly degenerate 
energy levels. 

In QuTiP there are two solvers that calculate the time 
evolution according to Eq. (|9]): odesolve numerically in- 
tegrates the set of coupled ordinary differential equations 
(ODEs), and essolve which employs full diagonalization. 
The odesolve and essolve solvers both take the same 
set of input parameters (as exemplified in Sec. |4]) and can 
easily be substituted for each other in a QuTiP program. 
For a quantum system with N states, the number of ele- 
ments in the density matrix is N'^ , and solving the master 
equation by numerical integration or diagonalization in- 
volves of use of superoperators of size x N"^. In the 
sparse matrix format, not all of the N"^ elements need to 
be stored in the memory. However, the time required to 
evolve a quantum system according to the master equa- 
tion still increases rapidly as a function of the system size. 
Consequently, the master equation solvers are practical 
only for relatively small systems: N < 1000, depending on 
the details of the problem. In Fig. |3] we show the scaling 
of the elapsed time for a typical simulation, here chosen to 
be the Heisenberg spin-chain 



M M-1 



(10) 



+ Jy<Jy<Jy 



Jz <^z <^z \ 



as a function of the size of the Hilbert space, for the two 
master equation solvers, as well as for two realizations of 
the Monte-Carlo solver mcsolve described the following 
section. In general, the exact time required to evolve a sys- 
tem depends on the details of the problem, but the scaling 
with system size is rather generic. The Monte-Carlo solver 
has superior scaling properties compared to the master- 
equation solvers, but due to the overhead from stochas- 
tic averaging, it is only for systems with a Hilbert space 
dimension around ~ 1000 that the Monte-Carlo solvers 
outperform the master equation. 
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Figure 3: ( Col or) The time required to evolve the Heisenberg spin- 
chain, Eq. | |lO| , as a function of the system size 2^^ where M is the 
number of spins, using the master equation ODE solver odesolve 
(blue), diagonalization via essolve (green), and the Monte-Carlo 
solver mcsolve with 250 (red) and 500 (cyan) trajectories, respec- 
tively. The dashed lines give the estimated calculation times extrap- 
olated from the data when the simulation could no longer fit in the 
computers memory (odesolve) , or the calculation became intractable 
(essolve). Here, the spin parameters are assumed to be identical 



with h = 2tt and Jx 



J y 



J. 



0.1 X 27r. Likewise, each spin 
has a dephasing rate given by 7 = 0.01. The initial state is given by 
= |l)l|0)2 • ■ ■ |0) M- Calculations were performed on a 2.8 Ghz 
quad-core computer with 24 Gb of memory. 



conditioned on the increase in information gained about 
the state of the system via the environmental measure- 
ments [8|. In general, this evolution is governed by the 
Schrodinger equation ([s]) with a non-Hermitian effective 
Hamiltonian 



cff 



ih 



2 E 



(11) 



where again, the C„ are collapse operators, each cor- 
responding to a separate irreversible process with rate 
7„. Here, the strictly negative non-Hermitian portion of 



Eq. (Ill gives rise to a reduction in the norm of the wave 



function, that to first-order in a small time is given by 

{%}){t -f 5t)\'4){t + 6t)) = 1-Sp where 



(12) 



3.2. Monte- Carlo trajectories 

Where as the density matrix formalism describes the en- 
semble average over many identical realizations of a quan- 
tum system, the Monte-Carlo (MC), or quantum-jump ap- 
proach [4] to wave function evolution, allows for simulating 
an individual realization of the system dynamics. Here, 
the environment is continuously monitored, resulting in 
a series of quantum jumps in the system wave function. 



and St is such that (5p <C 1. With a probability of remain- 
ing in the state \'ip{t + St)) given by 1 — Sp, the correspond- 



ing quantum jump probability is thus Eq. ( 12 1. If the en- 
vironmental measurements register a quantum jump, say 
via the emission of a photon into the environment [35 ] , or 
a change in the spin of a quantum dot |36j , the wave func- 
tion undergoes a jump into a state defined by projecting 
\ip(t)) using the collapse operator C„ corresponding to the 
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^1/2 



measurement 

mt+st)) = c,,\m)/{m\ciCn\my" ■ (13) 

If more than a single collapse operator is present in 



1.0 -I 



Eq ( 11 1, the probability of collapse due to the ith-operator 



Ci is given by 



(14) 



Evaluating the MC evolution to first-order in time is 
quite tedious. Instead, QuTiP uses the following algorithm 
to simulate a single realization of a quantum system [371 
ISHllSn]- Starting from a pure state |V'(0)): 

I: Choose a random number r between zero and one, 
representing the probability that a quantum jump oc- 
curs. 

II: Integrate the Schrodinger equation ([s]), using the 
effective Hamiltonian ( 11 ) until a time t such that the 
norm of the wave function satisfies {iP{t) {ipiT)) — 
at which point a jump occurs. 

Ill: The resultant jump projects the system at time r 



into one of the renormalized states given by Eq. ( 13 ). 



The corresponding collapse operator C„ is chosen such 
that n is the smallest integer satisfying 



> r 



(15) 



i=l 



where the individual P„ are given by Eq. (14). Note 



that the left hand side of Eq. (15) is, by definition, 
normalized to unity. 

IV: Using the renormalized state from step III as the 
new initial condition at time t, draw a new random 
number, and repeat the above procedure until the fi- 
nal simulation time is reached. 

3.2.1. Example: Single-photon cavity decay 

As an illustrative example, let us consider the evolution 
of a single-photon cavity Fock state in a non-zero ther- 
mal environment [101. The evolution of the wave func- 
tion |'0(O)) = |1) is governed by the effective Hamiltonian 



H, 



1 



off 



'-^iKa'^a- 



(16) 



with cavity frequency Wc, cavity decay rate k, and where 
is the steady state thermal occupation number. 



/th 



While the first term in Eq. ( 11 ) is responsible for the stan- 



dard unitary evolution of the cavity mode, the second and 
third terms give rise to random quantum jumps to lower 
and higher cavity photon numbers, respectively. When 
a jump occurs, the wave function of the system is pro- 
jected into a state corresponding to the collapse operator 



§ 0.5 

0.0 
1,0 

§ 0.5 

0.0 
1.0 

§ 0.5 

0.0 
1.0 

g 0.5 



0.0 





0.0 



0.2 0.4 
Time (sec) 



0.6 



Figure 4: (Color) (a) Monte-Carlo simulation showing the number 
operator expectation value for a single trajectory (blue) in the de- 
cay of a single-photon Fock state from a cavity coupled to a thermal 
environment, as demonstrated experimentally in Ref. 1401 . Here, the 
average photon number for the thermal environment is n = 0.063 
(black), while the decay rate k = 1/Tc is given by the cavity ring- 
down time Tc = 0.129. (b-d) Averages of 5, 15, and 904 trajecto- 
ries showing ensemble averaging toward the master equation solution 
(dashed-red) . 



Ci — y^(T^j-~(?Tyj^JjKa, yielding a decrease in the cavity oc- 
cupation number, or C2 = ("')th ■• which results in an 
increase. Here, the relative ratio of jumps corresponding 
to an increase in the cavity occupation number to those 
for decay is determined by the magnitude of A sin- 

gle realization of this evolution, showing a lone quantum 
jump, is presented in Fig. [4^. 

In addition to single quantum systems, when averaged 
over a sufficiently large number of identical system realiza- 
tions, the MC method leads to the same evolution equa- 
tion as the Lindblad master equation (|9| for a pure state 
density matrix ^ [8] . Therefore, the MC method may be 
used in any situation where the Lindblad master equation 
is valid, as discussed in Sec. |3.1[ However, for large quan- 
tum systems with Hilbert space dimension N ^ 1, the 
MC method is vastly more efficient than simulating the full 
density matrix given that only N elements are required to 
simulate a wave function, as opposed to the N'^ elements 
necessary in the ME approach. Although multiple trajec- 
tories are required, convergence to the ME result scales as 
■m~^ [311, where m is the number of trajectories simulated. 
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Figure 5: (Color) Scaling of the sum of the absolute deviation per 
time step, averaged over 10 simulations, for expectation values calcu- 
lated with the Monte-Carlo solver and compared to the ME solution, 
as a function of the number of trajectories. The error measure is cal- 
culated by integrating the absolute errors over the entire evolution 
time, and normalized by the number of time steps. The acceptable 
level of error varies from case to case, but in general the error measure 
should be much smaller than unity. A 1/m fit to the data (dashed) 
showing the predicted convergence rate is also presented. 




Figure 6: (Color) The time-evolution of a two-qubit system described 
by the Hamiltonian Eq. | |17[ ), which at time t/T = 1 and ideal con- 
ditions (dashed lines) transforms the qubit states in accordance with 
the i-SWAP gate. With qubit relaxation and dephasing (solid lines), 
the end-result deviates from the ideal i-SWAP gate. For these par- 
ticular parameters {g = 2it, Fi = 0.75, r2 = 0.5), the fidelity of the 
dissipative gate is 91%. 



In typical situations, between 250 and 500 trajectories are 
sufficient for errors smaller than a few percent, see Fig. [5] 
In Fig. [4] we show the convergence of the MC simulation 
to the ME solution for the single-photon cavity example 
as the number of trajectories averaged over is increased. 

4. Numerical calculations 

In this section we illustrate, via a number of examples, 
how quantum dynamical calculations are carried out using 
the QuTiP framework. The typical workflow for perform- 
ing a simulation with qutip is: 

I: Define the parameters that characterize the system 
and environment (if applicable). 

II: Create Qobj class instances representing the 
Hamiltonian and initial state of the system. 

Ill: For dissipative systems, define the collapse oper- 
ators as Qobj objects. 

IV: Evolve the system with a choice of evolution algo- 
rithm and output (e.g., operator expectation values). 

V: Post-process and visualize the data. 

Using the quantum object described in Sec. [2] and the 
solvers for the time-evolution of quantum systems de- 
scribed in Sec.[3j we can explore a diverse set of problems. 
The examples presented here are selected because they il- 
lustrate different features of QuTiP with a minimum of 
complexity. 



4-1. Fidelity of a two-qubit gate subject to noise 

To introduce how the evolution of a dynamical quantum 
system is calculated using QuTiP, let us consider a simple 
system comprised of two qubits that, during a time T = 
7r/4(7, are subject to the coupling Hamiltonian 

H ^ g {cr^ (g) a-a; + ay (g) ay) , (17) 

where g is the coupling strength. Under ideal conditions 
this coupling realizes the i-SWAP gate between the two 
qubit states [HJ |35] . This can readily be seen by evolving 
any initial state for the time T, and comparing the final 
state with the corresponding i-SWAP transformed initial 
state. We shall assume that the qubits are coupled with 
their surrounding environments, resulting in qubit energy 
relaxation and dephasing. 

Following the workflow outlined in the previous section, 
the QuTiP code for this problem can be organized in the 
following manner. First, define the numerical constants in 
the problem. For brevity, the code for this step has been 
omitted. Next, the Qobj instances for the Hamiltonian 
and the initial state may be defined as 

H = g * (tensor (sigmaxO , sigmaxO) + 
tensor (sigmayO , slgmayO)) 

psiO = tensor (basis (2 , 1) , basis(2,0)) 

To model qubit relaxation and dephasing, we define a list 
of collapse operators that later will be passed on to the 
ODE solver. For each qubit we append its associated col- 
lapse operator to this list (here called c_op_list) 

sml = tensor (sigmamO , qeye(2)) 

szl = tensor (sigmazO , qeye(2)) 

c_op_list . append(sqrt (gl * (l+nth) ) * sml) 

c_op_list . append(sqrt (gl * nth) * sml.dagO) 

c_op_list . append(sqrt (g2) * szl) 
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where the parameter nth is the number of 
environmentany-induced thermal excitations in the 
steady state. The collapse operators containing (T_ and 
al_ describe the qubit relaxation and excitation with 
the rates gl * (1+nth) and gl * nth, respectively, and 
the collapse operator models qubit dephasing. These 
lines of codes are repeated for the second qubit, with the 
appropriate change in the definition of the operators (i.e., 
the arguments in the tensor function are switched). 

At this point we are ready to let QuTiP calculate the 
time-evolution of the system. In the following example we 
use the master equation ODE solver odesolve. In addition 
to the Hamiltonian, initial state, and the list of collapse 
operator, we pass a list tlist to the solver that contains 
the times at which we wish to evaluate the density matrix. 

tlist = linspace(0, T, 100) 

rho_llst = odesolve(H, psiO, tlist, c_op_list, [] ) 
rho_final = rho_list [-1] 

If the last parameter is empty, as in this example, all 
QuTiP time-evolution solvers return the full density ma- 
trix (or state vector) corresponding to the times in tlist. 
Alternatively, a list of operators may be passed as last 
argument to the solver, in which case it will return the 
corresponding expectation values. 

Given the output of density matrices, we may now cal- 
culate the corresponding expectation values for selected 
quantum operators. For example, to calculate the excita- 
tion probability of the two qubits as a function of time, we 
may use the QuTiP function expect 

nl = expect(siiil .dagO * sml, rho_list) 
n2 = expect (siii2 . dagO * sm2, rho_list) 

Here, nl and n2 are now real NumPy arrays of expectation 
values, suitable for plotting or saving to file. 

Finally, to quantify the difference between the lossy i- 
SWAP gate and its ideal counterpart, we calculate the 
fidelity fidelity 

U = (-Ij * H * pi / (4*g) ) . expmO 
psi_ideal = U * psiO 

rho_ideal = psi_ideal * psi_ideal . dagO 
f = fidelity (rho_ideal, rho_f inal) . 

The results are shown in Fig. [6] where the expectation 
values for the two qubits, as a function of time, is plotted 
both with and without dissipation. The full code for this 
example is listed in [Appendix B.3 



1.00 



4-2. J aynes- Gumming model 

The same method used in the previous section can calcu- 
late the dynamics of the Jaynes-Cumming model Eq. ([2| . 
Only the definitions of the Hamiltonian, initial state, and 
collapse operators need to be changed to solve this prob- 
lem. For the Jaynes-Cumming model, the Hamiltonian 
and a possible initial state were given in Sec. [21 and we 
need only to define collapse operators before tne system 
can be evolved using one of the QuTiP solvers. The cavity 
and the atom relaxation rates are k and P, respectively. In 
this example, only the cavity is coupled to an environment 
with Boltzmann occupation number nth- We can write the 
collapse operators for the cavity 

a = tensor (destroy (N) , qeye(2)) 

c_ops . appendCsqrt (kappa * (l+n_th) ) * a) 

c_ops . append(sqrt (kappa * n_th) * a.dagO) 




Time 



Figure 7: (Color) Evolution of the Jaynes-Cummings Hamiltonian 
|[2| in a thermal environment characterized by {n)^^ = 0.75. Initially, 
only the atom is excited, but the atom-cavity coupling results in a 
coherent energy transfer between the two systems, a phenomenon 
known as vacuum Rabi oscillations. Here, the atom and the cavity 
are resonant, ujo = e = 2n, the coupling strength g/uio = 0.05, 
and the atom and cavity relaxation rates are •y / (u]o/2it) = 0.05 and 
K/(aJo/27r) = 0.005, respectively. 



and for the atom 

sm = tensor (qeye(N) , destroy(2)) 
c_ops . append (sqrt (gamma) * sm) 



Instead of having the solver return the state, as in Sec. 4.1 
we can request that the expectation values for a list of op 
erators be directly calculated at each time step. In this 
Jaynes-Cumming problem we are interested in the excita- 
tion number of the cavity and the atom, and as such can 
then define a list of expectation value operators 

expt_ops = [a.dagO * a, sm.dagO * sm] 

which may be passed as last argument to, for example, the 
odesolve solver. 



tlist = linspace(0, 10, 100) 

expt_list = odesolve(H, psiO, tlist, c_ops. 



expt_ops) 



The solver then returns a NumPy array expt_list of ex- 
pectation values. The result of this calculation is shown in 
Fig. [T] and the complete code is listed in [Appendix B.4[ 

4-3. Trilinear Hamiltonian 

To demonstrate the QuTiP Monte-Carlo (MC) solver, 
we consider the trilinear Hamiltonian that, in the interac- 
tion frame, may be written as 



H = ihK (afo^c^ - a) he) 



(18) 



consisting of three harmonic oscillator modes convention- 
ally labeled pump (a), signal ib) and idler (c) respectively, 
with the frequency relation, uja = uJi, + Wc and coupling 
constant K. This Hamiltonian is the full quantum gener- 
alization of the parametric amplifier |43| describing several 
quantum optics processes, including frequency conversion 
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[33] , the interaction of two- level atoms with a single mode 
resonant EM field [31] , and the modeling of Hawking ra- 
diation from a quantized black hole 46] . Here we suppose 
the pump mode is initially in a coherent state, while the 
signal and idler modes are in the ground state 



>a|0>JO), 



(19) 



As a system comprised of three harmonic modes, this 
model readily lends itself to MC simulation since the 
Hilbert space dimensionality increases exponentially with 
the number of initial excitations in the system (A^(0))^ = 
|a|^. For example, to accurately model an initial pump 
mode coherent state with \a\^ = 10 requires ~ 17 states, 
suggesting a minimum Hilbert space dimensionality of 
17'^ = 4913 for simulating Eq. (18 1; a value five times 



larger than what can typically be efficiently calculated us- 
ing the odesolve or eseries solvers (see Fig. [s]). 

In QuTiP, the Hamiltonian ( 18 ) may be expressed as 
{K=l) 

H=l j * (a*b . dag *c . dag () -a . dag () *b*c) , 

with the destruction operators for the pump, signal, and 
idler modes, a, b and c respectively, created via the tensor 
product 

a=tensor (destroy(N) ,qeye(N) ,qeye(N)) 
b=tensor(qeye(N) , destroy (N) ,qeye(N)) 
c=tensor (qeye (N) , qeye (N) .destroy (N) ) . 

In addition, we may define number operator expectation 
values and collapse operators in the same manner as pre- 
vious examples 

numO.numl ,nuiii2= [aO . dagO *aO ,al . dagO *al , a2 . dagO *a2] 

CO , CI , 02= [sqrt (2 . 0*gO) *aO , sqrt (2 . 0*gl) *al , sqrt (2 . 0*g2) *a2] 

mcsolve takes the same input arguments as odesolve, 
save for an additional argument necessary to specify the 
number of MC trajectories to simulate 

ntraj=500 

In Fig. M we plot the expectation values for the three 
modes oithe trilinear Hamiltonian (18 1, with correspond- 
ing damping rates, -fa — 0^ — tTS, 7c = 0.2, for the 
initial state given by Eq. ( 19 1 with a ~ \/lO 



psiO=tensor (coherent (N , sqrt (10) ) .basis (N . 0) .basis (N . 0) ) 
avgs=iiicsolve (H.psiO ,tlist ,ntraj . [CO. CI .02] . [numO .numl ,num2] ) 

Had we not defined any collapse operators, the evolution 
calculated by mcsolve reduces to the Schrodinger equation 
([5]). This evolution is also presented i n Fig, ([g]). Th e 
underlying QuTiP code may be found in [Appendix B.^] 

4.4- Landau-Zener transitions 

Landau-Zener transitions [47] are an interesting prob- 
lem that involves a quantum two-level system with a time- 
dependent energy splitting. The Hamiltonian is 

(20) 



H{t) = 



where A is the tunneling rate, and v is the rate of change in 
the bare qubit energy splitting. The Landau-Zener tran- 
sition theory analytically describes how the final state at 




Time 



Figure 8: (Color) Occupation numbers for the three modes of the 
triUnear Hamilto nian \18\ , averaged over 1000 trajectories, for an 
initial state Eq. | |l9| with a = \/lO. In this simulation, the envi- 
ronment is assumed to be at zero temperature. The pump (blue), 
signal (green), and idler (red) mode damping rates are 7a = 0.2, 
76 = 0.8, and 7c = 0.2, respectively. The closed system evolution 
(dashed-colors) is also presented, where the idler mode is omitted as 
its evolution is identical to that of the signal. 



i — >■ cx) is related to the initial state at t — >■ —00. In par- 
ticular, the probability of an adiabatic transition from |1) 
to |0) is given by the Landau-Zener formula 



P = 1 - exp - 



(21) 



Using QuTiP, we can easily integrate the system dynamics 
numerically, and obtain the state of the system for any 

intermediate value of t. 

The Landau-Zener problem differs from the previous 
examples in that the Hamiltonian is explicitly time- 
dependent. The QuTiP solvers odesolve and mcsolve 
both support time-dependent Hamiltonians. In order to 
specify an arbitrary time-dependent Hamiltonian, a call- 
back function may be passed as the first argument to the 
time-evolution solvers (in place of the Qobj instance that 
normally represents the Hamiltonian). The callback func- 
tion is expected to return the value of the Hamiltonian at 
the given point in time t, which is the first argument to 
the callback function. 

def hamiltoiiian_t (t , args) : 
HO = args[0] 
HI = args[l] 
return HO + t * HI 

In addition, a second argument passed to the callback 
function is a user-defined list of parameters. For perfor- 
mance reasons, it is appropriate to let this list contain 
pre-calculated Qobj instances for the constant parts of the 
Hamiltonian. For the Landau-Zener problem this corre- 
sponds to 

HO = delta/2.0 * sigmaxO 
HI = v/2.0 * sigmazO 

H_args = (HO. HI) 

The list of arguments for the Hamiltonian callback func- 
tion is then passed on to the time-evolution solver (as the 
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Figure 9: (Color) The occupation probability of the |1) (red curve) 
and |0) (blue curve) states of a quantum two-level system throughout 
a Landau-Zener transition. The solid black line is the final state 
according to the Landau-Zener formula Eq. | |21| l. The parameters 
used in this calculation are A = 0.5 X 2-k and v = 2.0 x 27r. 
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Figure 10: (Color) Bloch sphere representation of the Landau-Zener 
transition presented in Fig. [9] Here, the color of the data points 
corresponds to the evolution time. 



very last argument), along with the callback function itself 
(as first argument). 

expt_list = odesolve (haiiiiltoiiiaii_t , psiO, tlist, 

c_op_list, expt_list, H_args) 

The result of this calculation is shown in Fig. (|9]), which 
shows the intermediate dynamics of the system in terms 
of the occupation probabilities of the |0) and |1) states. 
The full code is shown in Appendix B.6[ Adding the 



operators sigmax, sigmay and sigmaz to expt_list, this 
evolution can also be visualized on the Bloch sphere using 
QuTiP's built in Bloch class, as demonstrated in Fig. [TO 



The QuTiP code for this figure is given in [Appendix B .7 



Although simple, this example illustrates the support for 
time-dependent Hamiltonians in QuTiP. It is a straightfor- 
ward exercise to implement an arbitrary time-dependence 
by changing the definition of the Hamiltonian callback 
function, or by modifying the time-independent part of the 
Hamiltonian to correspond to a more complicated quan- 
tum system. 

5. Performance 

As with any scientific simulation, the performance of 
the underlying numerical routines in QuTiP is an impor- 
tant factor to consider when choosing which software to 
implement in the analysis of the problem at hand. In sim- 
ulating quantum dynamics on a classical computer, this 
is especially important given that creating composite sys- 
tems using the tensor product leads to an exponential in- 
crease in the total Hilbert space dimensionality. Thus, it 
is beneficial to compare the performance of QuTiP to the 
other currently available quantum simulation packages. In 
this section we compare the simulation times of the master 



equation and Monte-Carlo solvers in QuTiP, to the those 
of the qotoolbox, as a function of Hilbert space size. 

To compare the QuTiP master equation solver, 
odesolve, to the qotoolbox equivalent (also called 
odesolve), we evaluate the coupled oscillator equation 

{n = i) 



(22) 



with uja = ^ 27r and coab — 0.1 x 2tt. In addi- 
tion, we consider the situation in which one resonator is 
damped with a corresponding dissipation rate g — 0.05. 
Here, the initial state of the system is the tensor prod- 
uct of Fock states |?A(0)) = \N)a\N - l)b, where N is 
the number of states in the truncated Hilbert space for 
each oscillator. The total time needed to simulate the 
dynamics over the time range t € [0, 10] is shown in 



Fig. (11). We see that the QuTiP solver easily outper- 



forms the qotoolbox as the Hilbert space dimensionality 
D = increases. For the largest dimensionality con- 



sidered in Fig. (11) D = 200, QuTiP is 4 times faster 



than the qotoolbox single-precision C-code implementa- 
tion, even though QuTiP is performing double-precision 
calculations with a small Python overhead. 



The trilinear Hamiltonian model from Sec. 4.3 provides 
a useful demonstration of the multiprocessing routines 
used by the QuTiP mcsolve function. Here, the indepen- 
dent MC trajectories are run in parallel, with the num- 
ber of simultaneous trajectories determined by the num- 
ber of processing cores. For the large Hilbert spaces asso- 
ciated with the trilinear Hamiltonian, the increased over- 
head generated from multiprocessing is overcome by the 
gains in running Monte-Carlo trajectories in parallel. In 



Fig. (12), we highlight these performance gains in simu- 
lating Eq. by plotting the computation time over a 
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Figure 11: (Color) Comparis on o f computation times, averaged over 
three trials, for solving Eq. l |22[ l in both QuTiP (blue) and the qo- 
toolbox (red) as a function of Hilbert space dimension. The shaded 
region highlights the increasing performance benefit from using the 
QuTiP solver as the dimensionality increases. Simulations were per- 
formed on a quad-core 2.8 Ghz processor. 



range of Hilbert space dimensions. For comparison, we 
also plot the times required for identical simulations us- 
ing the qotoolbox, which is limited to a single processor. 
For this simulation, the qotoolbox outperforms our QuTiP 
implementation for system sizes D < 500, even when mul- 
tiple processors are utilized. This is due to the Python 
overhead needed to implement the Monte-Carlo algorithm 
However, as shown in Fig. ([3|, 



discussed in Sec. 3.2 However, as shown m J^ig. us- 
ing the master equation is more appropriate for systems 
of this size. As with the master equation solver. Fig. ( |Tl| ), 
the benefits of using the QuTiP Monte-Carlo solver be- 
come appreciable as the system size increases. Even for a 
single-processor, the QuTiP mcsolve routine outperforms 
the qotoolbox after D w 1500, where the Python over- 
head is no longer a limiting factor. When using multiple 
processing cores, the Monte-Carlo solver performance gain 
nearly equals the number of processors available, the ideal 
situation. 



6. Conclusion 

We have presented a new, open-source framework for 
the numerical simulation of open quantum systems im- 
plemented in the Python programming language. This 
framework is suitable for a wide range of computational 
problems in quantum systems, including unitary and dis- 
sipative time-evolution, spectral and steady-state proper- 
ties, as well as advanced visualization techniques. In this 
work we have described the basic structure of the frame- 
work, the Qobj class, and the primary evolution solvers, 
odesolve and mcsolve. In addition, we have highlighted 
a number of examples intended to give the reader a flavor 
of the types of problems for which QuTiP is suitable. For 



Figure 12: (Color) Comparison of computation times, averaged over 
three runs, between QuTiP and the qotoolbox (dashed) in simulating 
the trilinear Hamiltonian from Eq.JTsl for an initial pump mode 
coherent state with expectation value (N)a = v^/4 as a function 
of the Hilbert space dimensionality D. The average performance 
enhancement from using multiple processors, as compared to the 
single-processor performance, is given in the legend. Ideally, this 
value should be equal to the number of processors. Computations 
were performed on a quad-core 2.8 Ghz processor. 



more in-depth documentation, and more elaborate exam- 



ples using the functions listed in Table A.2 we refer the 
reader to the QuTiP website [23]. There, one may down- 
load the latest version of this framework, as well as find 
installation instructions for the most common computer 
platforms. The version of the framework described in this 
paper is QuTiP 1.1.3. 

As with any initial software release, the performance of 
QuTiP can likely be improved in the future, with addi- 
tional portions of the code being optimized with Cython 
[15] or implemented in PyOpenCL [49] . 
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States 

basis / fock 
coherent 
cohereiit_dm 
f ock_dm 
qstate 
thermal _dm 

Operators 

qeye 

create 

destroy 

displace 

squeez 

mun 

sigmax 

sigmay 

sigmaz 

sigmap 

sigmam 

jmat 

Functions on states 

entropy_vii 

expect 

fidelity 

ket2dm 

liouvillian 

orbital 

ptrace 

qf unc 

simdiag 

tensor 

tidy up 

tracedist 

wigner 

Evolution 

correlation_es 

correlation_mc 

correlation_ode 

correlation_ss_es 

essolve 

mcsolve 

Odeoptions 

odesolve 

propagator 

pr opagat or _st eady st ate 
steadystate 

Utilities 

about 

Bloch 

clebsch 

demos 

eseries 

ode2es 

par for 

Qobj 

sphereplot 



Creates single basis (Fock) state in Hilbert space. 
Single-mode coherent state with complex amplitude a. 
Coherent state density matrix with complex amplitude a. 
Density matrix representation of a single basis (Fock) state. 

Tensor product state for any number of qubits in either the ground or excited states. 
Thermal state density matrix. 



Identity operator. 
Bosonic creation operator. 
Bosonic annihilation operator. 
Single-mode displacement operator. 
Single-mode squeezing operator. 
Number operator. 
Pauli spin- 1/2 operator. 
Pauli spin-1/2 ay operator. 
Pauli spin-1/2 (t^ operator. 
(T+ operator. 
(T_ operator. 
Spin-j operator. 



Von-Neumann entropy of a density matrix. 
Calculates the expectation value of an operator. 
Calculates the fidelity between two density matrices. 
Converts a ket vector to a density matrix. 

Assembles the Liouvillian super-operator from a Hamiltonian and a list of collapse operators. 

Calculates an angular wave function on a sphere. 

Partial trace of composite quantum object. 

Husimi-Q function of a given state vector or density matrix. 

Simultaneous diagonalization of commuting Hermitian operators. 

Calculates tensor product from list of input operators or states. 

Removes small elements from a quantum object. 

Trace distance between two density matrices. 

Wigner function of a given state vector or density matrix. 



Two-time correlation function using exponential series. 

Two-time correlation function using Monte-Carlo method. 

Two-time correlation function using ODE solver. 

Two-time correlation function using quantum regression theorem. 

State or density matrix evolution using exponential series expansion of ODE. 

Stochastic Monte-Carlo wave function solver. 

Options class for ODE integrators used by mcsolve and odesolve. 
ODE solver for density matrix evolution. 

Calculates the propagator U (t) for a density matrix or wave function. 
Steady state for successive applications of the propagator U{t). 
Calculates the steady state for the supplied Hamiltonian. 



Information on installed version of QuTiP and its dependencies. 
Class for plotting vectors and data points on the Bloch sphere. 
Calculates a Clebsch-Gordon coefficient. 
Runs built-in demos scripts. 

Exponential series representation of a time-dependent quantum object. 
Exponential series describing the time evolution of an initial state. 
Parallel execution of a for-loop over a single-variable. 
Class for creating user-defined quantum objects. 
Plots an array of values on a sphere. 



Table A. 2: List of user-accessible functions in QuTiP. Additional information about each function may be obtained by calling ^ functiori-namel' 
from the Python command-line, or by going to the QuTiP website 1231 . 
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Appendix B. QuTiP codes 

In this section we display the QuTiP codes underlying 
the calculations performed in generating Figures [2| |4] [6j 
[7j [8j |9j and [TOj For brevity, the code segments associ- 
ated with plotting have been omitted. The codes, in their 
entirety, may be viewed at the QuTiP website [25] , 



Appendix B.l. 



Figure [![■ 
Model 



Non-RWA Jaynes-Curamings 



from qutip import * 

## set up the calculation ## 

wc=1.0*2*pi# cavity frequency 

wa =1.0*2* pi # atom frequency 

N = 20 # number of cavity states 

g = linspaceCO, 2.5, 50)*2*pi # coupling strength vector 

## create operators ## 

a = tensor (destroy(N) , qeye(2)) 

sm = tensor(qeye(N) , destroy(2)) 

nc = a.dagO * a 

na = sm.dagO * sm 

## initialize output arrays ## 

na_expt = zeros (len(g) ) 

nc_expt = zeros (len(g) ) 

## run calculation ## 

for k in range (len(g) ) : 

## recalculate the hamiltonian for each value of g ## 

H = wc*nc+wa*na+g [k] * (a. dag()+a) * (sm+sm. dagO ) 

## find the groundstate ## 

ekets, evals = H.eigenstatesO 

psi_gnd = ekets [0] 

## expectation values ## 

na_expt [k] = expect(na, psi_gnd) # qubit occupation 
nc_expt[k] = expect (nc, psi_gnd) # cavity occupation 

## Calculate Wigner function for coupling g=2.5 ## 

rho_cavity = ptrace(psi_gnd,0) # trace out qubit 

xvec = linspace (-7 . 5 ,7 . 5 , 200) 

## Wigner function ## 

W = wigner (rho_cavity, xvec, xvec) 



Appendix B.2. Figure ^ Monte- Carlo relaxation in a 
thermal environment 

from qutip import * 

N=5 # number of basis states to consider 

a=destroy(N) # cavity destruction operator 

H=a.dag()*a # harmonic oscillator Hamiltonian 

psiO=basis (N, 1) # initial Fock state with one photon 

kappa=l . 0/0 . 129 # coupling to heat bath 

nth= 0.063 # temperature with <n>=0.063 

## collapse operators ## 

c_op_list = [] 

## decay operator ## 

c_op_list . append(sqrt (kappa * (1 + nth)) * a) 
## excitation operator ## 

c_op_list . append(sqrt (kappa * nth) * a.dagO) 
## run simulation ## 
ntraj=904 # number of MC trajectories 
tlist=linspace (0 ,0.6, 100) 

mc = mcsolve (H,psiO , tlist ,ntraj , c_op_list , [] ) 
me = odesolve (H,psiO , tlist , c_op_list , [a.dag()*a]) 
## expectation values ## 
exl=expect (num(N) ,mc [0] ) 

ex5=sum( [expect (num(N) ,mc [k] ) for k in range(5)] ,0)/5 
exl5=sum( [expect (num(N) ,mc [k] ) for k in range (15) ], 0) /15 
ex904=sum( [expect (num(N) ,mc[k]) for k in range (904) ] ,0)/904 



Appendix B.3. Figure^ Dissipative i-SWAP gate 

from qutip import * 

g =1. 0*2* pi # coupling strength 
gl = 0.75 # relaxation rate 

g2 = 0.05 # dephasing rate 

n_th =0.75 # bath temperature 

T = pi/(4*g) 

H = g * (tensor (sigmaxO , sigmaxO) + 

tensor (sigmayO , sigmayO)) 
psiO = tensor (basis (2 , 1) , basis(2,0)) 
c_op_list = [] 

## qubit 1 collapse operators ## 

sml = tensor (sigmamO , qeye(2)) 

szl = tensor (sigmazO , qeye(2)) 

c_op_list . append(sqrt (gl * (l+n_th)) * sml) 

c_op_list . append(sqrt (gl * n_th) * sml.dagO) 

c_op_list . append(sqrt (g2) * szl) 

## qubit 2 collapse operators ## 

sm2 = tensor(qeye(2) , sigmamO) 

sz2 = tensor(qeye(2) , sigmazO) 

c_op_list . append(sqrt (gl * (l+n_th)) * sm2) 

c_op_list . append(sqrt (gl * n_th) * sm2.dag()) 

c_op_list . append(sqrt (g2) * sz2) 

## evolve the system ## 

tlist = linspace(0, T, 100) 

rho_list = odesolve(H, psiO, tlist, c_op_list, [] ) 

rho_final = rho_list[-l] 

## calculate expectation values ## 

nl = expect (sml . dagO * sml, rho_list) 

n2 = expect (sm2 . dag * sm2, rho_list) 

## calculate the fidelity ## 

U = (-Ij * H * pi / (4*g)) .expmO 

psi_ideal = U * psiO 

rho_ideal = psi_ideal * psi_ideal . dagO 
f = f idelity(rho_ideal, rho_final) 



Appendix B.4- 



Figure ^ 
model 



Dissipative Jaynes- Cumming 



from qutip import * 

N = 5 # number of cavity states 

omegaO = epsilon = 2 * pi # frequencies 

g=0.05*2*pi # coupling strength 

kappa = 0.005 # cavity relaxation rate 

gamma =0.05 # atom relaxation rate 

n_th =0.75 # bath temperature 

## Hamiltonian and initial state ## 

a = tensor (destroy(N) , qeye(2)) 

sm = tensor (qeye(N) , destroy(2)) 

sz = tensor (qeye(N) , sigmazO) 

H = omegaO * a.dagO * a + 0.5 * epsilon * sz 

+ g * (a.dagO * sm + a * sm.dagO) 
psiO = tensor(f ock(N,0) , fock(2,l)) # excited atom 
## Collapse operators ## 
c_ops = [] 

c_ops . append (sqrt (kappa * (l+n_th)) * a) 

c_ops . append (sqrt (kappa * n_th) * a.dagO) 

c_ops . append (sqrt (gamma) * sm) 

## Operator list for expectation values ## 

expt_ops = [a.dagO * a, sm.dagO * sm] 

## Evolution of the system ## 

tlist = linspace(0, 10, 100) 

expt_data = odesolve (H, psiO, tlist, c_ops, expt_ops) 

Appendix B.5. Figure^ Trilinear Hamiltonian 
from qutip import * 

N=17 # number of states for each mode 
## damping rates ## 
g0=g2=0 . 1 
gl=0.4 
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alpha=sqrt(10) # initial coherent state alpha 
tlist=linspace(0,4,201) # list of times 
ntraj=1000#number of trajectories 
## lowering operators ## 
aO=tensor (destroy (N) ,qeye(N) ,qeye(N)) 
al=tensor(qeye(N) ,destroy(N) ,qeye(N)) 
a2=tensor (qeye(N) ,qeye(N) ,destroy(N)) 
## number operators ## 

n0,nl,n2= [aO . dagO *aO , al . dagO *al , a2 .dagO *a2] 
## dissipative operators ## 

CO , CI , C2= [sqrt (2 . 0*gO) *aO , sqrt (2 . 0*gl) *al , sqrt (2 . 0*g2) *a2] 
## initial state ## 

psiO=tensor(coherent(N, alpha) ,basis(N,0) ,basis(N,0)) 

## trilinear Hamiltonian ## 

H=lj*(aO*al .dag()*a2.dag()-a0.dag()*al*a2) 

## run Monte-Carlo ## 

avgs=mcsolve(H,psiO,tlist,ntraj , [C0,C1,C2] , [n0,nl,n2] ) 
## run Schrodinger ## 

reals=mcsolve (H.psiO , tlist , 1 , [] , [nO ,nl ,n2] ) 

Appendix B.6. Figure^' Landau-Zener transitions 

from qutip import * 

## callback function for time-dependence ## 
def hamiltonian_t (t , args) : 

HO = args[0] 

HI = args[l] 

return HO + t * HI 
delta = 0.5 * 2 * pi 
v=2.0*2*pi# sweep rate 
## arguments for Hamiltonian ## 
HO = delta/2.0 * sigmaxO 
HI = v/2.0 * sigmazO 
H_args = (HO, HI) 
psiO = basis(2,0) 
## expectation operators ## 
sm = destroy(2) 

sx=sigmax() ; sy=sigmay () ; sz=sigmaz() 
expt_op_list = [sm.dagO * sm,sx,sy,sz] 
## evolve the system ## 
tlist = linspace(-10.0, 10.0, 1500) 
expt_list = odesolve (haiiiiltonian_t , psiO, tlist, 
[] , expt_op_list , H_args) 

Appendix B. 7. Figure Block sphere representation of 
Landau-Zener transition 
Following the code from Appendix B.6 



import matplotlib as mpl 

from matplotlib import cm 

## create Bloch sphere instance ## 

b=Bloch() 

## normalize colors to times in tlist ## 
nrm=mpl . colors . Normalize (-2 , 10) 
color s=cm. j et (nrm (tlist) ) 

## add data points from expectation values ## 

b . add_points ( [p_ex [1] ,p_ex [2] , -p_ex [3] ] , 'm' ) 

## customize sphere properties ## 

b .point _color=list (colors) 

b . po int _mar ker = [ ' o ' ] 

b . point_size= [8] 

b.view=[-9,ll] 

b.zlpos=[l.l,-1.2] 

b.zlabel=['$\left |0\\right>_{f}$' , ' $\lef t I l\\right>_{f }$ ' ] 
## plot sphere ## 
b . showO 
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